# -*- coding: utf-8 -*-
#!/usr/bin/env python

import pylab, atpy, numpy, scipy, math

#my_params = atpy.Table('/home/jeffrey/Dropbox/isooctane/data/my_params_no_hfs.vot')
tbl_overall = atpy.Table('/home/jeffrey/Dropbox/isooctane/data/vizier_votable.vot')
#params_found = open('distance.out','w')
params_found_V = open('distance_V.out','w')
#params_curve = open('curve.out','w')
	
#mag1 and mag2 are the colours that you are interested in for the isochrone
#cut gives the magnitude you care about. This allows limiting to just the GB.
def dartmouth_isochrone_loader(filename,mag1,mag2,cut,distance_modulus,extinction):
	isochrones = open(filename,'r')
	col1=numpy.array([]);
	col2=numpy.array([]);
	for line in isochrones:
		if line[0]=='#':
			if line[0:4]=='#EEP':
				double=line.split()
				for i in range(len(double)):
					if double[i]==mag1:
						col1_index=i
						continue
					if double[i]==mag2:
						col2_index=i
						continue
			continue
		double=line.split()
		if float(double[col2_index])+distance_modulus+extinction<cut:
			col1=numpy.append(col1,float(double[col1_index]))
			col2=numpy.append(col2,float(double[col2_index]))
	return col1,col2

[col1,col2]=dartmouth_isochrone_loader('/home/jeffrey/Dropbox/isochrones/dartmouth/FeH161_121Gyr_BVRI.txt','B','V',16.5,13.7,0.3)
polycoeffs = scipy.polyfit(col1-col2+0.12,col2+13.7+0.3,5)
print polycoeffs
V_fit = scipy.polyval(polycoeffs, col1-col2+0.12)

LEID=tbl_overall.data['LEID']
B=tbl_overall.data['Bmag']
BmV=tbl_overall.data['B-V']

Vmag=B-BmV

params_found_V.write('LEID,BmV,V,distance\n')
#This does perpendicular distance from the isochrone
#for i in range(len(Vmag)):
	#distance_current=99
	#distance_best=99
	#distance_new=0
	#if math.isnan(BmV[i])==True:
		#continue
	#for BmV_curve in range(6000,16000):
		#BmV_curve=float(BmV_curve)/10000
		#distance_current=(math.sqrt(math.pow((BmV_curve-BmV[i]),2)+math.pow((scipy.polyval(polycoeffs, BmV_curve)-Vmag[i]),2)))
		#if abs(distance_best)>abs(distance_current):
			#distance_best=distance_current
			#if BmV_curve>BmV[i]:
				#distance_best=distance_best*-1
	#params_found.write(str(LEID[i])+','+str(BmV[i])+','+str(Vmag[i])+','+str(distance_best)+'\n')
	
#for BmV_curve in range(60,160):
	#BmV_curve=float(BmV_curve)/100
	#V_curve=scipy.polyval(polycoeffs, BmV_curve)
	#params_curve.write(str(BmV_curve)+','+str(V_curve)+'\n')
	
#This version just does the distance in (B-V) from the isochrone
polycoeffs_V = scipy.polyfit(col2+13.7+0.396,col1-col2+0.12,5)
for i in range(len(Vmag)):
	if math.isnan(BmV[i])==True:
		continue
	distance=BmV[i]-scipy.polyval(polycoeffs_V, Vmag[i])
	params_found_V.write(str(LEID[i])+','+str(BmV[i])+','+str(Vmag[i])+','+str(distance)+'\n')

pylab.figure(1)
pylab.plot(BmV, B-BmV, 'k*')
pylab.plot(col1-col2+0.12, col2+13.7+0.396, 'r-')
pylab.ylim([17,11])
pylab.xlim([-0.3,1.9])
#pylab.plot(col1-col2+0.12, V_fit, 'b-')
#pylab.savefig('isochrone.png')
pylab.show()



